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Abstract 

We propose a general parameter-free model identification tech- 
nique for a broad class of problems characterized by oscillatory dy- 
namics with a stable limit cycle using measurement data. The model is 
cast in the form of an autonomous descriptor system with an evolution 
equation for the dominant oscillation and with manifolds for the low- 
and high-frequency components. The descriptor system comprises the 
Landau equation, the mean-field model for a Hopf bifurcation, and 
more general Galerkin models of fluid flow as special cases. We develop 
and validate a variational data assimilation approach which allows us 
to identify the system by making assumptions only on the smoothness 
of the propagator. The proposed model identification technique is il- 
lustrated using transient vortex shedding in a wake flow as an example 
problem. It is demonstrated that this approach can be used to system- 
atically refine existing models, so that they describe more accurately 
available data. The article is written for practitioners working in the 
area of applied dynamical systems as the intended audience. 
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1 Introduction 

In this study, we propose a parameter-free model identification technique 
for oscillatory dynamics which broadens the spectrum of methods available 
for reduced-order modelling. Many real-life physical systems have high- 
dimensional state spaces reflecting the interaction of many subsystems or 
due to their spatially extended properties necessitating description in terms 
of systems of partial differential equations (PDEs). Yet, the behavior of their 
solutions may often be described using low- dimensional stable fixed points 
or stable limit cycles. Then, their dynamics may be approximated within 
a low-dimensional state space M N using evolution equation = f(a) for 
the state vector a = [ai, . . . , <in} t G Mr with propagator / : Mr —> M N , 
/ = [fu ■ ■ ■ > I'n] t - 

In some cases the propagator may be derived from the description of the 
full plant via first principles. In general, the reduced-order model is inferred 
from experimental or numerical data. Such reduced models are of paramount 
importance as test-beds for the development of our understanding of system 
dynamics. They are also useful as low-cost surrogates guiding optimization 
and real-time control design for expensive full-scale models. In all cases, typ- 
ical model identification is generally performed in three steps: (1) selection of 
the state space M N which is large enough to capture the behaviour of interest 
and at the same time sufficiently small to allow one to exploit the analyti- 
cal/numerical advantages of the surrogate plant; (2) structure identification 
of the propagator /, e.g., determination of the polynomial degree of its com- 
ponents fi, i — 1, . . . , N; and (3) parameter identification, e.g., inference of 
the polynomial coefficients from time- resolved trajectories t i— > a(t). The 
challenge of this approach is to find the right balance between the robust- 
ness of the model identification, requiring only a few tunable parameters, 
and a good accuracy for which a larger number of parameters is typically 
needed. Identification problems are often solved using variational techniques 
and this is the approach we will pursue in the present study. Such meth- 
ods have been developed for a broad range of problems in both the finite 
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and infinite-dimensional setting, including flow control in fluid mechanics 
P, data assimilation in meteorology (2j |3] and geophysics jl] to mention just 
a few application areas. The related problem of state estimation is usually 
solved using various filtering approaches such as the Kalman filter j5]. In 
this study, we show how the concepts concerning optimal reconstruction of 
constitutive relations developed in [51 [7] can be adapted to parameter-free 
propagator identification in a system characterized by oscillatory dynamics. 

While the idea of the proposed approach may be applicable to a broader 
class of problems, to fix attention, in the present investigation we focus 
on identification of a three-dimensional (3D) model with the state variable 



a(£) = [ai(t), a 2 (t), a 3 (t)] T e M 3 , where t > is time. Denoting 



:d 



where % is the imaginary unit, r := yaf+~a| is the amplitude of the fluc- 
tuations (":=" means "equal to by definition") and 9 := arctan(a2/ai) the 
corresponding phase, we make the following assumptions about the structure 
of the system 

Assumptions 1 (i). In the plane (01,02) the system exhibits an unstable 
fixed point at the origin and an attracting limit cycle, 



(ii). the dynamics is phase-invariant, i.e., the vector field depends only 



on 



(Hi), the state variable a% is "slaved" to a\ and 02, i.e., 03 = 03(01,02). 

As regards the time t, we will assume that t e [0, T] for some T > 0. As 
a general form of a dynamical system consistent with Assumptions [1] we will 
consider 



or, equivalent ly 

7ft 



at(t) 
a 2 (t) 

a 3 (t) 



f(t)= 9l {r{t))r(t), 
e\t)=g 2 (r(t)), 
as(*) = 9s(r(t)), 



( gi (r)I + g 2 {r)3) 
9M*)), 



at(t) 




~h 


_a 2 (t)_ 




A 



(2a) 

(2b) 
(2c) 



(3a) 
(3b) 
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where I 



1 
1 



and J 



-1 

1 



Subsequently, we will also use the 



notation £ := [a ly a 2 ] T and X := [0, r max ], where r max : = sup te[0)T] r(t). In 
([2]) and ([3]) the functions gi : X — > R, i — 1, 2, 3, are assumed sufficiently 
regular to make the systems well posed (the question of the regularity of 
functions gi, i = 1,2, 3, will play an important role in our approach and will 
be addressed in detail further below). Without loss of generality, we will also 
assume that 9 > 0. System fl3a|) . or fT2a| - (12bj) . is complemented with the 
initial condition, respectively, £(0) = £° := [a^, a^] 7 and r(0) = ||£°||, #(0) = 
ar ct an (a^/a?). Dynamical systems of the type ([2]) or ([3]) arise commonly 
as a result of various rigorous and empirical model-reduction strategies in 
diverse application areas such as fluid mechanics [5], thermodynamics [H] and 
phase transitions [10]. The main contribution of this work is development 
of a computational technique allowing one to reconstruct functions g±, g 2 
and #3 in ([2]) in a very general form based on some measurements. The 
key idea is to formulate a least-squares minimization problem in which one 
of the functions gi, i = 1,2, 3, is the control variable. Then, a variational 
gradient-based approach can be employed to find the optimal solution in a 
suitable function space. To fix attention, but without the loss of generality, 
in the present investigation we will consider system (T5]) as a reduced-order 
model of hydrodynamic instabilities in open shear flows obtained using a 
suitable Galerkin projection. Details of this problem are introduced in the 
next Section. In the following Section we develop our model identification 
approach, whereas in Section H] we present a number of computational results 
concerning identification of the model for the system considered in Section [2j 
Then, in Section we analyze the computational performance of our method 
and compare the results with the predictions of some standard models applied 
to the problem in question. Summary and conclusions are deferred to Section 
[6j whereas in Appendices |Aj [B] and [O we collect some technical results. 



2 Example Problem — Model Identification 
For a Vortex Shedding Instability 

In this Section we define a model identification problem associated with the 
transient two-dimensional (2D) cylinder wake which is a well-studied hy- 
drodynamic instability [TH [12], [13]. This flow is a representative example 
of phenomena characterized by the Hopf bifurcation with an unstable fixed 
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point and a stable limit cycle. Additional examples in this category include 
the Rossiter modes of the flow over a cavity [T3j and other shear flows [l_5~j. 
First, in Section |2~T] we describe the initial-boundary value problem for the 
infinite-dimensional Navier-Stokes system which is a system of coupled PDEs 
representing the conservation of mass and momentum in the motion of vis- 
cous incompressible fluid. Then, in Section 12. 2\ a low-dimensional Galerkin 
expansion is recalled which reduces the kinematic description down to three 
modes whose amplitudes serve as the state variables a\, 02 and a%. It will 
be demonstrated that the system governing these variables is in fact in the 
form 03]). Special forms of this system arising as models in numerous appli- 
cations are discussed in Section 12. 3\ whereas in Section 12.41 we describe the 
measurements used as the basis for the reconstructions. 

2.1 Cylinder Wake Flow 

The 2D flow around a circular cylinder is described here in the Cartesian 
coordinate system x := (x,y). The origin coincides with the center of the 
cylinder, the x-coordinate points in streamwise direction, while y represents 
the transverse coordinate. The velocity field u := (u,v) has components u 
and v aligned with the x-axis and y-axis, respectively. The pressure field is 
denoted p. The Newtonian fluid considered here is characterized by uniform 
density p and kinematic viscosity v. The flow properties depend on the 
Reynolds number Re := U Djv. Here, we consider the case of Re = 100 
which is far above the critical Reynolds number of 47 characterizing the 
onset of vortex shedding (i.e., the Hopf bifurcation) [HI [16], and far below 
the transitional Reynolds number of 187 which marks the onset of three- 
dimensional instabilities |17j . 

In the following, all quantities are assumed to be non-dimensionalized 
with the oncoming flow velocity U, the cylinder diameter D and density p. 
The flow is considered in a rectangular domain surrounding the cylinder 



n := {(x, y) : x 2 + y 2 > 1/4 A -5 < x < 15 A \y\ < 5} 



(4) 



and its evolution is described by the Navier-Stokes equation 



V • u = 

d t u + u ■ Vw = — Ait — Vp 

Re 



in f2, 



(5a) 



in f2. 



(5b) 
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The boundary conditions for the velocity comprise the no-slip condition at 
the cylinder, the free-stream condition at the inlet boundary, the free-slip 
condition at the lateral boundaries and the no-stress condition at the down- 
stream outflow boundary. The initial condition for the velocity at t = is 
the solution u s of the steady Navier-Stokes system (i.e., the unstable fixed 
point) perturbed by the real part of the corresponding most unstable stability 
eigenmode u\(x) 

u(x,0) = u 8 {x) + 0.02u£(se). (6) 

The field u*(x) is normalized to have unit L2^l) norm. 

The initial-boundary- value problem is numerically integrated with a finite- 
element method based on an unstructured grid and details are provided in 
[T8] . Figure [27T1 depicts three flow snapshots corresponding to the initial con- 
dition, an intermediate transient state and the flow approaching the limit 
cycle. 

2.2 Galerkin Expansion 

Below we review approximate modelling approaches typically employed to ob- 
tain low-dimensional descriptions of the dynamics described by (jSJ), cf. [TU] . 
The quantities corresponding to the limit cycle will be denoted with the su- 
perscript V and in the present problem in which the limit cycle is stable 
we will therefore have r° = r max . The transient wake is usually character- 
ized by the superposition of a time-varying symmetric base flow u B with 
the length of the recirculating region decreasing with time as the vortex 
shedding develops and an antisymmetric oscillating field u' representing the 
vortex shedding. The change of the base flow is accurately captured by a 
single "shift mode" u A [20], while the oscillation field can be approximated 
by the first two Proper Orthogonal Decomposition (POD) modes Uj, i — 1, 2, 
computed at the limit cycle [211 UB] • Some information about the construc- 
tion of orthogonal bases using the POD approach is presented in Appendix 
lAl The resulting truncated Galerkin expansion thus takes the form 

u(x,t) w u B (x,t) + u'(x,t), (7a) 
u B (x,t) = u s (x) + a A (t) u A (x), (7b) 
u'(x,t) = ai(t) Ui(x) + a 2 (t) u 2 (x). (7c) 

One typically ignores deformations of the oscillatory modes during the tran- 
sient. These deformations do not significantly alter the model predictions 
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Figure 1: Flow snapshots at (a) t — (initial condition), (b) t = 24.8 
(intermediate transient state), and (c) t = 99.6 (periodic solution at the limit 
cycle). The domain shown represents the computational domain Q which also 
coincides with the region where the Galerkin expansion is defined. The flow 
patterns are visualized using streamlines (which are defined as the level sets 
of the streamfunction ip : Q — > M related to the velocity components via 

u = at and v = —jt-). The cylinder is indicated by the black circle. 

ov ox ' j j 
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and do not have any effect on the proposed model identification approach. 
On the other hand, inclusion of this effect in our model would significantly 
complicate the study and is outside the scope of this investigation. Iden- 
tifying a 3 = cxa and assuming phase-invariant behavior, we note that the 
evolution of the Galerkin expansion coefficients is governed by a system in 
the form Q, provided that gi(0) > 0, corresponding to an unstable fixed 
point at the origin, and g\{r°) = and dgi/dr\ r=r o < 0, corresponding to a 
locally attracting limit cycle at r = r° > 0. Moreover, from the assumptions 
made in Section [fl it follows that V r >o #2(7) > 0. Introducing a number 
of further simplifications one arrives at two well-known reduced models, the 
Mean-Field Model and the Landau Model. Since they provide a point of 
reference for our test problems, we describe them briefly below. 



The mean-field model of oscillatory flow instabilities [22] explicates the mech- 
anism of amplitude saturation in the mean-field deformation. This deforma- 
tion is quantified by the amplitude a a of the shift mode, cf. (I7bj) . which is 
slaved to the fluctuation level r. For the soft (supercritical) Hopf bifurcation, 
the resulting evolution equations read 



The ordinary differential equations (I8"a|) and (I8bj) correspond to the Navier- 
Stokes equation linearized around the time- varying base flow with as the 
order parameter. The third algebraic equation ( )8c|) can be thought of as the 
Reynolds equation linking the mean flow to the Reynolds stress generated 
by the fluctuations described by the first equation. 

Evidently, equations particular case of (EJ) with gi(r) = cr 1 — 

Pa o>a "c 2 ■> Qiij) — w i + 7a °a r2 an d <?3( r ) — «a r2 - The parameters a-y, Ui, 
Pa, 7a, O!^ of (jSJ) may be derived from the Galerkin approximation described 
in Section I2.2[ or identified from the solutions of the Navier-Stokes system 



2.3 Mean-Field Model 



f(t) 

e(t) 

ClA(t) 



[cri - /3 A a A (t)] r(t), 

^i + 7aoa(*), 
<r 2 (t), 



(8a) 
(8b) 
(8c) 
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via suitable fitting. The periodic solution of (jHJ) is given by 




(9b) 



(9a) 



oi //3a- 



(9c) 



At the limit cycle the mean flow is predicted to have a vanishing growth 
rate g%. This marginal stability property of mean flows has been conjectured 
by Malkus [23] and is corroborated by the global stability analysis of the 
Navier-Stokes system [21] . Derivation details and a further discussion of the 
properties of the mean-field model can be found in [T8j . 

The Landau equation for the supercritical Hopf bifurcation is a corollary 
to the mean-field model and constitutes a prototype evolution equation for 
self-amplified, amplitude-limited oscillations. It is obtained by substituting 
flHcD in (Iga MgEfr : 



Here, we require <7i,/3,Ui > to ensure a stable limit cycle with a positive 
angular velocity in the (01,02) plane, whereas 7 may vanish, be positive or 
negative. Evidently, equations (JTOj) correspond to (j8a]) - (j8bj) with (3 = «a/3a 
and 7 = a A 7A- 

The Landau equation (fTUj) may also be derived from the first principles, 
e.g., using a center-manifold reduction method, see, e.g., [23 ESI EZ] • The 
parameters may be identified from solutions of the Navier-Stokes system 
with 0i and oj\ obtained as the real and imaginary part of the most unstable 
eigenvalue associated with the Hopf bifurcation, whereas the nonlinearity 
parameters /3 and 7 can be inferred from the post-transient amplitude and 
frequency, r° and u°, respectively. 

2.4 Measurements 

The state in our approximate model ([2]) is characterized by three time- 
dependent mode amplitudes: a±, a<i and 03 and as "measurements" we will 
consider the functions a\(t), a-i (t) and a&(t) which are obtained for t G [0, T] 



r 







0i r — (3 r 3 , 
uji + 7 r 2 . 



(10a) 
(10b) 
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by solving initial-boundary- value problem (EJ)— (EJ) followed by the reductions 
described in relations (JTj). In the next Section we introduce a computational 
approach allowing one to optimally identity the functions gi{r), <?2( r ) and 
gs(r) in a suitable class, so that the predictions of system fl2]) best match 
in the least-squares sense the data di, a 2 and 5a • We note that, alterna- 
tively, a system of evolution equations for a 2 and a% may be obtained by 
substituting ansatz f[Taj) into momentum equation fl5b|) . and some compar- 
isons between such an approach and the results obtained using the model 
identification method developed here will be drawn in Section 15.31 

3 Computational Approach 

The task of identifying the functions gi, i = 1,2,3, such that the output 
of system fl2]) matches certain "measurements" is an example of an inverse 
problem j4]. What makes this problem somewhat different from typical in- 
verse problems is that the functions sought have the form of "constitutive 
relations", in the sense that they depend on the state variables (i.e., the de- 
pendent variables in the problem), rather than the independent variables. 
More specifically, in the problem considered here gi, i = 1,2,3, depend on 
r = \/ af + a\ as opposed to t. Non-parametric formulations of such inverse 
problems have received some attention in the context of systems described 
by PDEs [HI 128], but we are not aware of similar approaches applied to the 
phase-space description of dynamical systems. 

3.1 Formulation of Optimization Problem 

We will look for functions elements of the Sobolev space 

H 1 (X) of continuous functions with square- integrable gradients on X which is 
the "ident inability" region defined in Section [TJ i.e., the interval spanned by 
the state variable r(t) during the system evolution, see also [6]. Some addi- 
tional remarks concerning the regularity of functions gi and g^ are presented 
in Appendix [Bj In the reconstruction problems considered in this study, the 
boundary behavior of the reconstructed functions will have to be restricted in 
different ways. This reflects the fact that system ([2]) with the reconstructed 
functions gi, i = 1, 2, 3, should exhibit specific behavior in some distinguished 
regions of the phase space, namely, at the equilibrium r = and at the limit 
cycle r = r°. There are several possible choices and the boundary conditions 
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we adopt will reproduce the behavior described by the mean-field model in- 
troduced in Section 12.31 More specifically, 

• at the origin r = 

— the Jacobian of the right-hand side (RHS) in equation (12 aj) should 
be given by <?i(0) which is a priori unknown and is to be deter- 
mined as a part of the solution of the reconstruction problem; on 
the other hand, the Jacobian of the RHS in equation f|2bl) should 
vanish; this is achieved when the derivatives of both functions g± 
and g 2 are set to zero at the origin, i.e., 



dr 



0. 



1,2, 



(11) 



at the limit cycle r 



the RHS of 
i.e., 



2aj) should vanish resulting in the vanishing of gi, 



9i(r 



0. 



;i2) 



— in regard to equation ( I2bl) . we will prescribe a given slope G > 
of the RHS, i.e., 

^(r)U = G. (13) 

We note that, while the boundary behavior described by (JTT l) -(fr3~ l) at the 
equilibrium and the limit cycle is the same as in the mean-field model (cf. Sec- 
tion [23]), the behavior of the constitutive relations g\ and g 2 for intermediate 
values of the state magnitude < r < r° can be arbitrary and will be deter- 
mined using our optimal reconstruction procedure. No restrictions are placed 
on the boundary behavior of function g 3 , cf. ( )2c|) . 

A convenient way to solve inverse problems is to formulate them as suit- 
able optimization problems [21]. For each of the functions gi, i = 1, 2, 3, we 
thus define the corresponding cost functional Ji(gi) : H l {X) — > R as 



[r(t)-f(t)} dt, 



T 



D i9(t) 



dt, 



IdMt)) ~ «a(*)] 2 dt, 



(14a) 
(14b) 
(14c) 
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where r(t) := ^Jai(t) 2 + a 2 (t) 2 , 9{t) := arctan(a 2 (t)/ai(t)) and a&(t) are the 
"measurements" obtained as described in Section 12.41 In (j!4ap and (114bl) 
the functions r(t) and 6(f) are related to g\ and g 2 via system ([2]). The 
optimal reconstructions g~i, §2, and ^3 are defined as solutions of the following 
problems 

(PI) gt := argmin 9ieHl(I)j *_ gi{r)lr=0=0 , 9l(r c )=0 ^i(<7i), (15a) 
(P2) 32 := argmin 92eH1(I)> ^i^o, -^(OU^G ^(#2), (15b) 
(P3) # 3 := argmin g3eH i (I) Js^s)- (15c) 

Below we describe in detail a gradient-based approach to solution of Prob- 
lem (PI). Problem (P2) has a similar structure and the solution method is 
essentially the same with some modifications described hereafter. In Prob- 
lems PI and P2 it is assumed that the functions, respectively, #2 arid g\ are 
fixed. Problem (P3) has a different structure warranting a separate solution 
approach which will be described further below. 

3.2 Gradient-Based Approach to Solution of Problems 

(PI) and (P2) 

The minimizer gi is characterized by the first-order optimality condition [30] 
requiring the vanishing of the Gateaux differential 
Ji(9i,g[) ■= lim^oe -1 [Ji(gi + eg[) - Ji(gi)), i.e., 

Viewer), J^(r)| r=0 =o, 9i(r°)=o JKhiQi) = 0, (16) 

where g[ is an arbitrary perturbation direction. The (local) minimizer can 
be computed with the following iterative procedure 

^^f'-r^VJ!^), n = l,..., 

(i) ^ ' 

9i = 9i, 

where g® represents the initial guess, n denotes the iteration count and VJ\ : 
X — > R is the gradient of cost functional J\. The length of the step is 
determined by solving the following line minimization problem 

r W = argmin T>0 J x [gf ] - rVjrf)) (18) 
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which can be done efficiently using standard techniques such as Brent's 
method [HI]. For the sake of clarity, formulation (I17j) represents the steepest- 
descent method, however, in practice one typically uses more advanced min- 
imization techniques, such as the conjugate gradient method, or one of the 
quasi-Newton techniques [H2] • Evidently, the key element of minimization al- 
gorithm ( 1T7|) is the computation of the cost functional gradient Vj7i. It ought 
to be emphasized that, while the governing system ([3]) is finite-dimensional, 
the gradient Vj7i is a function of the state magnitude r and as such represents 
a continuous (infinite-dimensional) sensitivity of cost functional Ji{gi) to the 
perturbations g[ = g[(r). The fact that the control variable gi, and hence 
also the gradient Vj7i, are functions of the state variable, rather than the 
independent variable (Figure^), will result in cost functional gradients with 
structure rather different than encountered in typical optimization problems 
for differential equations (see [7] for some related questions arising in PDE 
optimization problems). 

In order to identify an expression for the gradient Vj7i, we proceed by 
computing the Gateaux differential of the cost functional 

J 1 (g 1 ;g[)= [r - r] r> \ 9l] g[) dt = — £ T £' dt = ( V J 1 (g 1 ), g[) 

Jo Jo r x ' x \ x ) 

(19) 

where we used the identity r' = £ T £'/ r and the last equality is a consequence 
of the Riesz representation theorem [33J with (-,-)x{x) denoting an inner 
product in the Hilbert space X(X) (to be specified later) of functions defined 
on X. The perturbation variable £' is a solution of the following perturbation 
problem (see Appendix [C] for a derivation) 

?(t) = [g 1 {r(t))I + I£(t) (VgMt))) T 

+g 2 (r(t))3 + J £{t) (Vg 2 (r(t))) T ] £'(t) + I£(t) g[ 
=: A($(t))£ + It(t)9v (20a) 
C'(0) = 0, (20b) 

r i T 

where = §^ , i = 1,2. We note that Gateaux differential f TT9]) 

is not yet in the form consistent with the Riesz representation, since the 
perturbation g[ does not appear in it as a factor, but is hidden on the RHS 
in perturbation equation (j2Uaj) . A standard technique to convert Gateaux 
differential f[T9l to the Riesz form is based on the adjoint variable : 
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Figure 2: (a) Schematic indicating the dependence of (thick red dashed- 
dotted line) the "constitutive relation" g i: i = 1,2, 3, on the state magnitude 
r with the plane (ai,a 2 ) representing the phase space; the thick blue solid 
line marks a sample trajectory (,2(t)] T , t e [0,T], of system (|3~aj) with 

the black circle representing the state at time to! ( D ) schematic illustrating 
the relation between the integration variables dt and dr, cf. ( 124")) ; the blue 
solid line represents the system trajectory C; the state £(t) (marked with a 
black point on the trajectory) together with the corresponding adjoint state 
£*(£) carry information necessary to evaluate cost functional gradient Vj7i(r), 
where r = \\£{t)\\, cf. ((26]) 



Optimal Model Identification 



15 



[0,T] -»■ M 2 . Taking the inner product (in R 2 ) of £*(t) with equation f l20al) . 
integrating over [0, T] and then integrating by parts we obtain 

= £ (D T {i' - [<?i(r)I + 1 1 (V^(r)) T + g 2 (r) J + J £ (V<? 2 (r)) T ] £' - £ dt 
= £ (O r { - r - [9i (r)I + U (V^(r)f + ^ 2 (r)J + J £ (V<? 2 (r)) T ] T £* } (ft 



T 

T 



(OS' - / (O j m)g'idt. 

t=0 



Defining the adjoint system as 



(21) 



-4* (t) = [Am)f V® + £, (22a) 

r 

C(T) = 0, (22b) 
we reduce relation ( l2Tj) to 

^i;^)= [ (e) T m)g'idt. (23) 

Jo 



We note that, although already appears as a factor in expression ( l2"3l . 
this expression is still not in the Riesz form, since the integration is with 
respect to the time dt, whereas in the inner product (v)at(z) defining the 
Riesz representer integration is with respect to the measure dr defined on the 
interval I (connection between the different integration variables is illustrated 
schematically in Figure [2b). The two variables are related via the following 
transformation 



£2 + £| ^ dr= ^ ^ - ^ ^ = dt, (24) 

where we used the identities d£i = f\dt and <i£ 2 = fzdt, cf. Denoting 
the trajectory in the state space C := {u te [ ,r] €(t) £ R 2 }, and combining 
(123]) with (124"]) we obtain the expression 

f f^*1 T I£ f r ' max (£*) T I£ 

Jl(gi;g[) = / , ~ f } —, \ gi(r)dr = / * 6{r)dr (25) 
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which is already in the required Riesz form. In (1251) the line integral over the 
contour C and the definite integral over the interval X are equal, because for 
dynamical system f l3a|) points on the contour C with the magnitude r G (0, r°) 
are unique, so that the map r — > | \£(t) \ = r) G C is one-to-one (in more 
general situations when this is not the case, or when the reconstructed func- 
tion depends on more than one state variable, e.g., both £i and £2 here, the 
change of variables needed to obtain the Riesz form will be more complicated 
and one has to employ more general techniques such as those developed in 

mm)- 

While this is not the gradient we will use in actual computations, we will 
first obtain an expression for the L 2 gradient which in the next subsection 
will be used as the basis for the calculation of gradients defined in the Sobolev 
space H l (T). Thus, setting X = £ 2 (Z) in {IS]), we obtain from (125]) 

Wi(r) = & V reI - (26) 

+ <S2j2 

Expression ( |26]) is validated computationally in Section 15.11 

As regards Problem P2, the optimality condition takes the form, cf. ( TT6|) 
and lfI5]) . 

V ^etfi(x), ^(r)| r=0 =o, A fli(r) | r=r0=G Ji(tosQ= sm(9-9)6'dt 

J 

sin(#-#)^— ^cft = 0, 
r 

(27) 

where we used the identity 6' = r 2 £ T J Following the same steps as 
described above, we obtain an expression for the cost functional gradient in 
the form (|2"6"1) . however, the adjoint system satisfied by has now a different 
source term on the RHS 

-t(t) = [Am)f C(t) + Sm{6 r ~ 9) J t(t), (28a) 
C(T) = 0. (28b) 



In regard to Problem P3, using change of variables (I24p . we can rewrite 
([He]) as 

1 /" r r , , , 2 



^3(^3) = > / T-rv n ^ ~ dr - ^ 
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Thus, in optimization problem (I15cjl we look for a function g% G H 1 ^) 
which is as close as possible (in a weighted L 2 topology) to a given function 
^3 £ -^(X), This problem, in fact, does not have a solution because of the 
density of the function space H 1 (X) in L 2 (X), cf. [31]. However, it is possible 
(and satisfactory from the application point of view) to "solve" problem (I15cl) 
approximately by finding a #3 G H 1 (X) such that Jz{gs) is sufficiently small. 
Such an approach is described in Section 13.31 

3.3 Sobolev Gradients 

In this Section we describe how Sobolev gradients V ffl Ji G H l (X), i = 
1,2, used in minimization algorithm ( 1TT1) for Problems PI and PI can be 
obtained from fl25|) . In addition to enforcing smoothness of the reconstructed 
functions, this formulation allows us to impose the desired behavior at the 
endpoints of the interval X, cf. ffTTl) - ffT3l . via suitable boundary conditions. 
We begin by defining the H 1 inner product on X as 

z lZ2 + e 2 -^-^-dr, (30) 

where I G 1R is a parameter with the meaning of a "length scale". It is 
well known [35] that extraction of cost functional gradients in the space 
H l with the inner product defined as in (|30|) can be regarded as low-pass 
filtering the L 2 gradients with the cut-off wavenumber given by t~ x . As 
regards the behavior of the gradients V H J at the endpoints of the interval 
X, we can require the vanishing of either the gradient itself or its derivative 
^:(V H J\ and the boundary conditions we prescribe correspond to relations 
(|TT]) - f|T3|) introduced as a part of the formulation of optimization problems 
fll5al) - fll5bl) . As regards the boundary data at r = r° (i.e., at the limit cycle), 
in Problem P2 we use ^:V H J 2 { r )\r=r° — which ensures that the property 
'dj92( r )\r=r° — G of the initial guess g 2 remains unchanged during iterations 
(fej). 

Identifying expression ( 125]) with inner product (l30j) . cf. ( TT9]) . integrating 
by parts and using the boundary conditions mentioned above we obtain the 
following elliptic boundary- value problem on X defining the Sobolev gradient 
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v H1 j 

( X " VH1 J = V ^ in (°. r °)> ( 31a ) 

^V H 'j = at r = 0, (31b) 

dr 

(PI) : V R1 J ) 

a , > = at r = r°, (31c) 

where the expression for V L2 J7" is given in (j26p . 

As concerns Problem P3, we propose to reconstruct G H 1 (X) directly 
(i.e., without iterations) by solving the following problem 



d 2 
dr 





= o 3 (r) 


in (0,r°), 


d t 








= 


at r = 0, 


9i 


= a 3 (r°) 


at r = r° 



(32a) 

(32b) 
(32c) 

which, expect for the boundary condition at r = r°, has an identical structure 
as ( IHTj) . The superscript in p| indicates dependence of the solution on the 
parameter i. We note that as £ — > the left-hand side (LHS) in f)32a|) 
approaches the identity transformation which means that ||^| — 03 ||i 2 (2:) 0, 
so that also Jz{gi) — > 0, as £ — > 0. Since solutions of system fl32|) are not 
defined for t = 0, we will obtain our approximate reconstruction as §3 := g| 
for some small value of I. Results concerning model identification for the 
system described in Section [2] are presented in Section HJ whereas in Sections 
15.11 and 15.21 we analyze certain computational aspects of the method. 



4 Results 

In this Section we present results concerning the solution of model identifi- 
cation problems PI, P2 and P3, cf. f ll5al) - fll5cl) . for the system introduced 
in Section [2j Motivated by practical considerations, we make the following 

Assumptions 2 In the solution of Problems PI and P2 we set in system 
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(GkD 

(PI) : 92 = 0, (33a) 

(P2) : 9l = 9l , (33b) 



which means that in the reconstruction of g<i we use the best available es- 
timate of g\ obtained from the solution of Problem PI. The choice of g?, 
has no effect on the solution of Problem PI and hence without loss of gen- 
erality we can adopt (I33al) . As the gradient descent algorithm in Prob- 
lems PI and P2 we use the Polak-Ribiere version of the nonlinear con- 
jugate gradient method [32] in which the "momentum" term is reset to 
zero every 20 iterations and iterations (1171) are declared converged when 

(Ji(9i n+1) ) ~ Ji(gl n) ))/Ji(9 { i n) ) < iO -7 , i = 1,2. In the system described in 
Section [2] the limit cycle is characterized by r° = 2.3, whereas the length 
of the time window is chosen as T = 70 which is long enough to allow the 
transient to settle on the limit cycle (see Figure [5|i below). We have suc- 
cessfully solved Problems PI, PI and P3 using different combinations of 
numerical parameters, and the parameters used to obtain the results pre- 
sented in this Section are summarized below. Systems ()3]) and (122"]) were 
solved using MATLAB subroutine ode45 with an adaptive time-stepping. 
Unless stated otherwise, the integrals defined on the interval [0, T], cf. (114)) . 
were discretized using = 500 equispaced points. The interval X was dis- 
cretized using N% = 75 equispaced points, and boundary- value problems ( l31~j) 
and ( 132]) were approximated using the second-order finite differences. The 
length-scale parameter appearing in (13T|) and (132|) was I = 1.0 in Problems 
PI and P2, and I = 0.1 in Problem P3. The initial condition £° for system 
(J2D and the initial guesses g® and g\ for reconstruction algorithm (|T7|) must 
be chosen so that the magnitudes t G [0, T], span the entire interval X, 

as otherwise the sensitivities (gradients) cannot be properly defined for all 
values of r (we refer the reader to [6] for a discussion how this limitation can 
be overcome in some cases). The best available initial guesses incorporate 
our prior knowledge about the system behavior at the origin and at the limit 
cycle within the structure of mean-field model (jHJ). We thus used 

g°(r) = 0.151 - 0.151 ' , g° 2 (r) = 0.886 + 0.15 (34) 
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respectively, for Problems PI and P2. It is clear that initial guesses (1341 
satisfy properties fj 1 X H — (1 X 3 H with G = 0.224. As the initial condition £° for 
([3]) we used a small perturbation around the fixed point at the origin. 

Our main results are presented in Figure [3] where we show the optimal 
reconstructions §i(r), i = 1,2,3, and compare them against the left-hand 
sides of equations (f2"a|) - (j2"cj) . all shown as functions of the state magnitude r. 
For completeness, the LHS of equations (|2ap — (]2b|) are shown as functions of 
time t G [0, T] in Figure 0] (in the case of gi, cf. Figures [3^ andHJ the LHS of 
equation (I2a| is additionally divided into r). In Figures [3K,b we also indicate 
the initial guesses (1341) which were used to obtain the optimal reconstructions 
via iterations (1T7I) . We reiterate that these initial guesses represent in fact 
the RHS of mean-field model ([HD with the best available choice of the param- 
eters. We see in Figures [3] that, as expected, the reconstructed constitutive 
relations gi(r), i = 1,2,3, smoothly approximate the left-hand sides of the 
corresponding equations evaluated using the measurements. Considered as 
functions of r, these left-hand sides are multi- valued which is a consequence 
of the fact that, due to the oscillations of the measurement data at the limit 
cycle (see Figure Eh below), the map t — > r(t) is not single- valued. Based 
on the reconstructions g± and cj2 we can obtain estimates of two important 
quantities, namely, the growth rate of the instability at the origin given by 
4- [gir] r=0 = gi(0) = 0.1576, and the oscillation frequency at the limit cycle 
given by §2^°) = 1.130. Finally, in Figure Owe compare the outputs from 
system ([2]) obtained using the mean-field model employed as initial guesses 
( )34|) and the reconstructions cji(r), i = 1,2, 3, against the corresponding mea- 
sured quantities (as regards the time-history of the state variables, we do not 
show ci2(t), as it has qualitatively very similar behavior to ai(t) already shown 
in Figure Eb). In Figures |5£i,b (see, in particular, the insets) we note that the 
evolution of r(t) and ai(t) obtained using the optimal reconstructions g± and 
g 2 is much closer to the measured quantities than the evolutions computed 
using the initial guesses g® and g® based on the mean-field model ([8]). We 
remark, however, that the measurements f(t) shown in Figure |5^l reveal some 
high-frequency oscillations which are not captured by the trajectory r(t) ob- 
tained using the optimal reconstruction §i. These oscillations reflect a phase 
dependence in the behavior of the solutions of the original Navier-Stokes sys- 
tem 05]), an effect which is by construction excluded from ansatz f l2a|) . As a 
consequence, g\ can convey phase-averaged information only. We will return 
to this problem again in Section 15.31 As regards the results shown in Fig- 
ure [5b, we see that, while the optimal reconstruction g^ir) is quite smooth 
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(a) 



IN 




0-9 !|J A 



(b) 




(c) 



Figure 3: (Blue solid lines) the optimal reconstructions of the constitutive 
relations (a) gi(r), (b) <?2( r ) and (b) g%{r)] (red dashed lines) the corre- 
sponding initial guesses (a) g1(r) and (b) ^(r); (black symbols) values of 
(a) r _1 (dr/dt)\ r ( ti ), (b) (d8 / dt)\ r ( ti ) and (c) aA\r(ti), computed based on the 
measurement data at the time instants t«, i = 1, . . . , N. 
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Figure 4: Dependence of the LHS in equations (I2a]) - (l2b]l . i.e., (red solid 
line) r~ l (dr / dt)\ r (t) and (blue dotted line) (d9 / dt)\ r (t) evaluated based on 
the measurement data on time t. 



(Figure Eb), the quantity g 3 (r(t)) exhibits oscillations absent in the original 
measurement data CI3 (t). This effect as well is a consequence of the lack of 
phase- dependence in ansatz ( |2a|) and ( l2cl ). 



5 Discussion of Computational and Physical 
Aspects 

In this Section we analyze a number of computational and physical modelling 
aspects of the proposed approach which can be important in applications. We 
begin by examining the accuracy of the cost functional gradients in Section 
15. H followed by a study of the robustness of iterations flTjl in Section 15.21 
and conclude with some insights about physical modelling in Section 15.31 

5.1 Validation of Gradients 

A key element of optimization algorithm (TlTl) are the cost functional gradi- 
ents and a standard approach to their validation consists in computing the 
directional Gateaux differential JKgug'i), i — 1,2, for some arbitrary per- 
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Figure 5: Comparison of (a) the state magnitude r, (b) state variable a 2 and 
(c) state variable a 3 as functions of time obtained from (red dashed line) the 
measurement data, (blue solid line) solution of system ((2]) using the optimal 
reconstructions g±, §2 and #3, and (black dotted line) solution of system 
02]) using initial guesses (134)) ; in Figures (a) and (b) insets are included to 
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turbations g[ in two different ways, namely, using a finite-difference approx- 
imation (with step size e) and using the inner product of the adjoint-based 
gradient with the perturbation g[, namely Riesz representation f|T9|) . and then 
examining the ratio of the two quantities, i.e., 

for a range of values of e. If the gradient V L2 j7i(r) is computed correctly, 
then for intermediate values of e, «j(e) will be close to the unity. Remarkably, 
this behavior can be observed in Figures [6^,b corresponding to Problems PI 
and P2 over a range of e spanning about 8 orders of magnitude. The quantity 
shown in Figures [6^,b is log |ftj(e) — 1|, i = 1, 2, which represents the number 
of significant digits to which the two ways to evaluate Jl{g%\ g'j) in (|35p agree. 
Furthermore, we also observe that refining the resolution Nt of the time 
interval [0, T] yields values of Kj(e) closer to the unity. The reason is that 
in the "optimize-then-discretize" paradigm adopted here such refinement of 
the discretization leads to a better approximation of the continuous gradient 
fl26l) . As can be expected, the quantities «j(e) deviate from the unity for 
very small values of e, which is due to the subtractive cancellation (round- 
off) errors in finite-differencing, and also for large values of e, which is due 
to the truncation errors, both of which are well-known effects. 



5.2 Computational Robustness of the Proposed Ap- 
proach 

In this Section we focus on the effect that the choice of initial guess g®, 
i — 1,2, has on the reconstructed functions g\ and g-i- We note that, given 
the nonlinearity of system (j2J), optimization problems PI and P2 may be 
nonconvex and optimality conditions (fT6l) and fl27|) characterize minimizers 
which are only local. Thus, different initial guesses may in principle give rise 
to different reconstructions and convergence to a global minimum cannot 
be a priori assured. We investigate this issue in Figures [Th. and [7b where 
we show the reconstructions obtained, respectively, in Problems PI and P2 
using different initial guesses. As regards Problem PI, we note in Figure [TJi 
that accurate reconstructions are obtained using even relatively poor initial 
guesses g®. On the other hand, in Figure [7b we see that in Problem P2 the 
reconstruction fails for a less accurate initial guess g®- These two examples 
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Figure 6: Diagnostic quantities (a) log |/«i(e) — 1| evaluated for Problem PI 
and (b) log|K 2 (e) — 1| evaluated for Problem P2, cf. ( I35p . as functions of 
log e obtained using different discretizations of the time interval with (blue 
circles) N T = 50, (red squares) N T = 500 and (black triangles) N T = 5000; 
in all cases the perturbation direction is g\ = — r 3 , i — 1,2. 



are representative of the behavior we generally observed in our calculations 
and we conclude that Problem PI appears more robust with respect to the 
choice of the initial guess than Problem P2. We also noted that in both 
problems convergence tends to be more sensitive to the values assumed by 
the initial guesses g\ and g\ at r = and r = r° than to their behavior 
for intermediate values of r. Results from Figure [7] are corroborated by 
the corresponding histories of the cost functionals in Figures [H^ and |Sb. 
We note that in the case of the poorest initial guess in Problem P2, cost 
functional Ji{g\ ) reveals hardly any decrease with the iterations at all. 
While in the cases of successful reconstructions the cost functionals Ji(g^) 
and ^(g^) drop over several orders of magnitude, they never attain values 
lower than 0(1O~ 3 ). This is a consequence of the phase-dependent behavior 
of the measurements which cannot be resolved using ansatz in the form fl2]), 
see the inset in Figure |5k, 

Inverse problems of the type considered here often tend to be ill-posed, 
in the sense that small perturbations to the data, for example due to noise, 
may result in significant changes in the computed solution. Suitable regular- 
ization, for instance, using Tikhonov's technique [U [29] , may be required to 
stabilize the solution procedure in such situations. To focus attention, in the 
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(a) (b) 

Figure 7: (Solid lines) reconstructions g^ i = 1,2, and (dashed lines) corre- 
sponding initial guesses g® in the solution of (a) Problem PI and (b) Problem 
P2. Reconstructions are marked with the same color as the corresponding 
initial guesses. Initial guesses and reconstructions discussed in Section @] are 
marked in black. 
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(a) (b) 

Figure 8: Decrease of the cost functionals with iterations in (a) Problem 
PI and (b) Problem P2 for the cases studied in Figure [7] (with the same 
color-coding). 
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present study we concentrated on the structure of the gradients and did not 
investigate the effect of noise on the reconstructions, hence such regulariza- 
tion was not necessary. We refer the reader to [6] for a thorough analysis of 
regularization applied to a related reconstruction problem. 

5.3 Physical Interpretation of the Results 

In this Section we propose some physical interpretation of the numerical 
reconstruction results from Section HI We note in Figures [3] and [5] that the 
identified phase-invariant oscillation model ([2]) is in fact in remarkably good 
agreement with the data obtained from the solution of the Navier-Stokes 
system (jHJ). A small difference between the model and the data is visible 
as wiggles due to the second harmonic present in the measurements which 
violates the phase-invariance assumed in our model ansatz. This difference 
can be easily removed by a simple pre-processing of the measurement data. 
Referring to Appendix |AJ we note that the POD eigenvalue Ai, representing 
the variance of a±, is larger than the eigenvalue A2 which represents the 
variance of a<i- The following rescaling transformation 



ensures equipartition of energy in the new variables a\ and a 2 while con- 
serving the total energy in both modes. This rescaling effectively removes 
the second harmonics in a^, % = 1,2, which could be used as new inputs for 
the reconstruction. However, we refrained from applying (|36|) in the com- 
putations reported in Section @] in order to show the power of the proposed 
identification method to deal with data which is not perfectly matched by 
the model. 

Another observation concerning Figure [3] is the significant deviation of 
reconstructed functions g\ and §2 from the parabolic mean-field relations 
OH]). Evidently, higher-order corrections, such as r 4 , r 6 , etc., are required 
for a better agreement between the identified propagators i = 1,2, and 
the expansions used in the mean-field model. The information which higher- 
order terms ought to be included in the model as opposed to an a priori fixed 
polynomial expansion used typically in model identification is therefore the 




(36b) 



(36a) 
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unique advantage of the proposed identification strategy. We note that odd 
powers of r can be excluded by phase-invariance considerations. 

These results also shed light on the validity of mean-field model OH]). 
Initially, the mean-field theory [361 EZ] was derived to be valid near the onset 
of the oscillation only. We probe the applicability of the model by applying 
it at a Reynolds number 100 which is more than twice the critical value of 47. 
Hence, the deviation of gi and §2 from f[TU|) does not invalidate the mean-field 
theory. One reason for this deviation is the change of the structure of the 
vortex street during the transient. The optimal oscillatory modes deform 
from the stability eigenmodes into the POD modes while the fluctuation 
center moves upstream and the frequency and wavenumber increase [T£| 13"%]. 
Similarly, the mean- field correction (the shift mode, cf. ( 17b |) ) changes during 
the transient [20] which has a noticeable effect on the mean- field model [39] . 

Finally, the results concerning the identified descriptor system are also 
relevant to the empirical 9-dimensional Galerkin model accounting for the 
base- flow variation and for the first four harmonics [18]. The initial expo- 
nential growth of the first harmonic is limited by the base-flow variation 
(which reduces the production of fluctuation energy) and, to a lesser extent, 
by the energy transfer from the first into higher harmonics. The energy 
transfer may be accounted for by an energy-dependent eddy viscosity in the 
mean- field system. Under certain assumptions (see, e.g. jlQl chap. 3]), a 
generalized Landau equation 

f = a"ir — (3r 3 — 7r 3 (37) 

can be derived, where <j\ denotes the growth rate near the fixed point r = 
and P, 7 characterize the damping from the 0-th and from higher harmonics, 
respectively. By comparing carefully the 9-dimensional Galerkin model with 
the identified phase-invariant system it may be therefore possible to deter- 
mine the values of /3 and 7, or even to correct the powers of the new terms 
in ( 13~Tj) . A complete derivation and an in-depth discussion of this problem is 
outside the scope of the present study. 

6 Conclusions and Future Directions 

A new method for model identification has been proposed and validated. 
We depart from the traditional approach of (1) characterizing the propaga- 
tor of the dynamical system in a parameter space and then (2) performing 
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a parameter identification. Thus, arbitrary polynomial expansions of the 
propagator may be performed a posteriori (following the solution of the op- 
timal reconstruction problem) at a practically vanishing cost. In addition, 
the performance of parametric models may easily be assessed and, if neces- 
sary, improved by introducing additional terms motivated by the form of the 
reconstructed constitutive relation. 

The method is applied to a three-dimensional descriptor system with 
three a priori undetermined relations describing the fluctuation growth, fre- 
quency and mean-field correction as functions of the fluctuation energy. As a 
benchmark problem we chose the onset of laminar von Karman vortex shed- 
ding behind a circular cylinder. Results of a direct numerical simulation are 
transcribed into the mode amplitudes of a minimal 3-state Galerkin model 
[T8] which are then captured with remarkable accuracy by our identified de- 
scriptor system. The form of the reconstructed system features a noticeable 
departure from the mean-field models which may be corrected by introducing 
in the latter higher-order terms. 

As regards future research directions, we remark that a surprisingly large 
set of modelling and control problems can be cast in a similar form of function 
identification of a descriptor system 



For reasons of simplicity, let us assume that function / is known and that 
function g needs to be determined. If b characterizes the slow modes, then 
b = g{a) represents the inertial manifold to be identified from a given sys- 
tem trajectory [27J. If b represents high-frequency components such as the 
parameters of a subgrid turbulence representation, e.g., the eddy viscosities 
[JT], then their functional dependence on the state variable a may also be 
inferred with our approach. On the other hand, if b denotes the actuation 
amplitudes, as in numerous wake stabilization studies [HJ H3j HU 05J 06], 
then b = g(a) represents a full-state feedback control law. This control law 
may as well be identified from desired trajectories t H- a(t). 

In summary, the proposed parameter-free modelling approach has unique 
advantages for a broad class of modelling and control problems. We are 
actively exploring these opportunities in the context of more complex fluid 
flow applications. 




(38) 
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A Proper Orthogonal Decomposition 

In this Appendix we describe the Proper Orthogonal Decomposition (POD) 
employed in Section [2721 to construct a low-dimensional model from the sim- 
ulation data. It is closely related to other techniques of data analysis known 
as the Principal Component Analysis or, in the discrete setting, the Singular- 
Value Decomposition. The starting point are snapshots of the velocity field 
u m (x):= u(x,t m ), m — 1, . . . , M, x e f2, cf. @. These snapshots are sam- 
pled at times t m uniformly spaced over one period of oscillation and form 
a statistically representative ensemble for the considered first and second 
moments. 

The goal is to construct a 'least-order' Galerkin expansion 

N 

U(X, t) = Uo(x) + ^ U i( X ) + U rcs(x, t) (39) 

i=l 

with base mode u$ and N space-dependent expansion modes (basis functions) 
Ui(x) with the corresponding mode amplitudes a;(t) which will result in a 
minimum- norm average residual it rcs of the snapshot ensemble (see, e.g., 
[471 140| HE]). The base mode is necessary so that Galerkin expansion ( )39l) 
satisfies inhomogeneous boundary conditions for arbitrary mode amplitudes. 
For instance, expansion ( 1391) captures the prescribed oncoming flow velocity 
regardless of the values of dj. 

The Galerkin expansion and its residual are embedded in the Hilbert 
space L 2 (f2) of square- integrable vector fields. The inner product of two 
elements v, w G L 2 (fl) is defined by 

( v ^ w )l 2 (u) : = J v-wdx, (40) 

Q 

where '•' denotes the standard Euclidean inner product and dx an infinites- 
imal volume element of the domain Q. The associated norm thus is 

||«lk(fi) := yJ(u,u) Lm - (41) 

We search for empirical modes Ui, i = 0, . . . , N, where N < M — 1, which 
minimize the average residual of the Galerkin expansion of the snapshots 

N 

u m = u + Y,<^ + < s , m = l,...,M, (42) 

j=i 
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with optimal mode amplitudes a™ 1 in the sense of the L2 norm, i.e 

M 



= J7rY^ IKsllL(n)= min - ( 43 ) 



i2 



m=l 

This problem is solved by the snapshot POD [49] and the base mode is the 
mean of the snapshots 

M 

-o:=^J> m (44) 

m=l 

The POD modes arise from the correlation matrix C := (C mn ) mn=1 M of 
the snapshot fluctuations 

C mn.= 1/ m_ n_ \ (45) 

M\ /L 2 (n) 

We note that C is a positive semi-definite Grammian matrix. Hence, the 
eigenvalue problem 

Ce^A^, z = l,...,M, (46) 

yields an orthonormal set of real eigenvectors = [ej, . . . , ef f ] , ef e 3 - = 8y, 
with non-negative eigenvalues which can be ordered as 

Ai > A 2 > ... > A M = 0. (47) 

The last equality arises from the fact that M vectors span a subspace of 
maximum dimension M — 1. Hence, M snapshots define only M — 1 POD 
modes and the corresponding amplitudes. These are given by 

1 M 

= -= ( W ™ - Uo ) , ar = ^A~Me7\ * = 1,...,M-1. (48) 

* * m=l 

The POD modes form an orthonormal basis in Z/2(fi), t*j)i a (n) = ^j? 
z,j = 1,...,N, while their amplitudes have vanishing means = 0, i = 
1, . . . , N, and diagonal second moments ai a] = A, 5^, i,j = 1, . . . , N. The 
average is to be understood in terms of the snapshot ensemble (see, e.g., 
(143]) ). The eigenvalue Aj can be interpreted as twice the fluctuation energy 
contained the i-th mode. 
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B Regularity of Reconstructed Function g\ 
versus Existence and Uniqueness of Solu- 
tions to Equation ( )2aJ ) 



By the one-dimensional embedding result H l {X) — C°'^(I) [51], we note 
that the reconstructed function will be Holder-continuous with A = 1/2. 
Thus, it will not meet the assumptions of the Picard-Lindelof theorem [50] . 
and in principle will ensure existence only, without uniqueness, of solutions 
of equation ( |2al) . In order to ensure the Lipschitz-continuity (A = 1) of 
gx, which would also guarantee uniqueness of solutions of ( |2a|) . we would 
need to reconstruct g\ as an element of Sobolev space H 2 (T), because then 
H 2 (T) — C 0,1 (X). While there are no fundamental difficulties here (we 
would need to replace inner product (130]) with the corresponding definition 
in H 2 (X)), we will refrain from this in the actual computations in Section H] in 
order to keep the approach as simple as possible. Nevertheless, in all problems 
we treated with the proposed approach the reconstructed functions possessed 
the Lipschitz regularity which was verified a posteriori by performing suitable 
grid-refinement studies. 



C Derivation of Perturbation Equation ( 120a 



In this Appendix we present a derivation of perturbation equation (I20ap . 
Obtaining this equation is made somewhat more involved by the fact that 
the perturbation variable g± is itself a function of the state magnitude, i.e., 
gi = gi(r). We assume here that only g\ is perturbed while g 2 remains 
fixed, with the opposite case leading to essentially the same calculations. 
By substituting, respectively, g\ = g\ a and gi = gu into equation (|3a| . 
we obtain £ a = (gi a (r a ) I + g 2 (r a ) J) £ a and = (fll&(r&) I + gi(r b ) J) £ 6 , 
where £ a := £(<?i a ) and := £(gib) are the corresponding solutions and 
r a '■= |£a|> r b '■= Taking the difference of these two equations and 
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defining |' := | a — |&- we obtain 

£' = 9la(r a ) I la + 92(r a ) J la - 01&(r&) I lb ~ S^b) J lb 

= 0i a (r o ) I |a + ^(^a) J |a + ^ib(r a ) I | a - 9ib(r a ) I la - fi'lb(rb) I lb - 02(rft) J lb 

= flf'(r a ) I la + #lb(r a ) I la - gib{n) I lb + ^(^a) J la - flbfo) J |b, 



A 



13 



(49) 

where we also set </(•) = <7i a (-) — fl'ib(-) in the first term on the RHS. As 
regards the terms denoted A in (149]) . they are transformed as follows using 
the fundamental theorem of calculus for line integrals and the change of 
variables |(s) = | 6 + s (| a - | 6 ) for s E [0, 1] 

A = g lb (r a )ie a -gib(r b )ie b = / V^d) |) d| 



jf V 4 F(| 6 + S |') ds) I', 



(50) 



where we also denoted F(|) := <?ib(|) |. The integrand expression on the 
RHS in (150]) is then expanded in the Taylor series around |' = 

■srftfa + sC) = ^(Ib) + or-aT^te) ^ s + (l^'l 2 ) ' ^ * = L 2 . 

(51) 

where the component notation was used for clarity with Fi, and |£ de- 
noting the components of vectors F, | and |'. Plugging expansion ( 15T1) into 
expression ( 150]) we obtain 



(52) 



A = V c (<7 16 (|) I) I' [ ds + 0(\$f) 
£=£b Jo 

= ~g lb (r b )I + U b (Vg lb (r b )) T ] |' + 0(||'| 2 ). 

Noting that term B in (1491) transforms in an analogous way to A in (l50l - (!52|) . 
using these results in (149]) . assuming smallness of |' and dropping terms of 
order quadratic and higher we finally arrive at perturbation equation ( I20a|) . 



References 



[1] M. D. Gunzburger. Perspectives in Flow Control and Optimization. 
SIAM, 2003. 



Optimal Model Identification 



35 



[2] I. M. Navon. Practical and theoretical aspects of adjoint parameter esti- 
mation and ident inability in meteorology and oceanography Dynamics 
of Atmosphere and Oceans, 27:55-79, 1997. 

[3] E. Kalnay. Atmospheric Modeling, Data Assimilation and Predictability. 
Cambridge University Press, 2003. 

[4] A. Tarantola. Inverse Problem Theory and Methods for Model Parameter 
Estimation. SIAM, 2005. 

[5] R. F. Stengel. Optimal Control and Estimation. Dover, 1994. 

[6] V. Bukshtynov, O. Volkov, and B. Protas. On optimal reconstruction 
of constitutive relations. Physica D, 240:1228-1244, 2011. 

[7] V. Bukshtynov and B. Protas. Optimal reconstruction of material prop- 
erties in complex multiphysics phenomena, (submitted), 2011. 

[8] J. Dusek, P. Le Gal, and Ph. Fraunie. A numerical and theoretical 
study of the first Hopf bifurcation in a cylinder wake. Journal of Fluid 
Mechanics, 264:59, 1994. 

[9] L. D. Landau and E. M. Lifshitz. Statistical Physics. Pergamon Press, 
1980. 

[10] I. R. Yukhnovskii. Phase Transitions of the Second Order — Collective 
Variables Method. World Scientific, 1987. 

[11] CP. Jackson. A finite-element study of the onset of vortex shedding in 
flow past variously shaped bodies. J. Fluid Mech., 182:23-45, 1987. 

[12] M. Morzyhski, K. Afanasiev, and F. Thiele. Solution of the eigenvalue 
problems resulting from global non-parallel flow stability analysis. Corn- 
put. Meth. Appl. Mech. Enrgrg., 169:161-176, 1999. 

[13] B. R. Noack and H. Eckelmann. A global stability analysis of the steady 
and periodic cylinder wake. J. Fluid Mech., 270:297-330, 1994. 

[14] C.W. Rowley and D.R. Williams. Dynamics and control of high- 
Reynolds number flows over open cavities. Ann. Rev. Fluid Mech., 
38:251-276, 2006. 



Optimal Model Identification 



36 



P. J. Schmid and D. S. Hennigson. Stability and Transition in Shear 
Flows. Springer- Verlag, New York, 2001. 

A. Zebib. Stability of viscous flow past a circular cylinder. J. Engr. 
Math., 21:155-165, 1987. 

H.-Q. Zhang, U. Fey, B. R. Noack, M. Konig, and H. Eckelmann. On 
the transition of the cylinder wake. Phys. Fluids, 7(4):779-795, 1995. 

B. R. Noack, K. Afanasiev, M. Morzyhski, G. Tadmor, and F. Thiele. A 
hierarchy of low- dimensional models for the transient and post-transient 
cylinder wake. J. Fluid Mech., 497:335-363, 2003. 

P. Holmes, J. L. Lumley, and G. Berkooz. Turbulence, Coherent Struc- 
tures, Dynamical Systems and Symmetry. Cambridge University Press, 
1996. 

G. Tadmor, O. Lehmann, B. R. Noack, and M. Morzyhski. Mean field 
representation of the natural and actuated cylinder wake. Phys. Fluids, 
22(3):034102-1..22, 2010. 

A. E. Deane, I. G. Kevrekidis, G. E. Karniadakis, and S. A. Orszag. 
Low-dimensional models for complex geometry flows: Application to 
grooved channels and circular cylinders. Phys. Fluids A, 3:2337-2354, 
1991. 

J.T. Stuart. Nonlinear stability theory. Ann. Rev. Fluid Mech., 3:347- 
370, 1971. 

W.V.R. Malkus. Outline of a theory of turbulent shear flow. J. Fluid 
Mech., 1:521-539, 1956. 

D. Barkley. Linear analysis of the cylinder wake mean flow. Europhysics 
Lett, 75:750-756, 2006. 

H. Haken. Synergetics, An Introduction. Nonequilibrium Phase Tran- 
sitions and Self- Organizations in Physics, Chemistry, and Biology. 
Springer, New York, 3rd edition, 1983. 

J. Guckenheimer and P. Holmes. Nonlinear Oscillations, Dynamical 
Systems, and Bifurcation of Vector Fields. Springer, New York, 1986. 



Optimal Model Identification 



37 



A. N. Gorban and I. V. Karlin. Invariant Manifolds for Physical and 
Chemical Kinetics. Number Vol. 660 in Lecture Notes in Physics. 
Springer- Verlag, Berlin, 2005. 

G. Chavent and P. Lemonnier. Identification de la non-linearite d'une 
equation parabolique quasilineaire. Applied Mathematics and Optimiza- 
tion, 1:121-162, 1974. 

C. R. Vogel. Computational Methods for Inverse Problems. SIAM, 2002. 

D. Luenberger. Optimization by Vector Space Methods. John Wiley and 
Sons, 1969. 

W. H. Press, B. P. Flanner, S. A. Teukolsky, and W. T. Vetterling. 
Numerical Recipes: the Art of Scientific Computations. Cambridge Uni- 
versity Press, 1986. 

J. Nocedal and S. Wright. Numerical Optimization. Springer, 2002. 

M. S. Berger. Nonlinearity and Functional Analysis. Academic Press, 
1977. 

R. A. Adams and J. F. Fournier. Sobolev Spaces. Elsevier, 2005. 

B. Protas, T. Bewley, and G. Hagen. A comprehensive framework for 
the regularization of adjoint analysis in multiscale pde systems. Journal 
of Computational Physics, 195:49-89, 2004. 

J.T. Stuart. On the non-linear mechanics of hydrodynamic stability. J. 
Fluid Mech., 4:1-21, 1958. 

J. Watson. On the non-linear mechanics of wave disturbances in stable 
and unstable parallel flows. Part 2. the development of a solution for 
plane Poiseuille flow and for plane Couette flow. J. Fluid Mech., 9:371- 
389, 1960. 

[38] G. Tadmor, O. Lehmann, B. R. Noack, L. Cordier, J. Delville, J. -P. Bon- 
net, and M. Morzyhski. Reduced order models for closed-loop wake con- 
trol. Philosophical Transactions of the Royal Society A, 369(1940):1513- 
1523, 2010. 



Optimal Model Identification 



38 



[39] M. Morzynski, W. Stankiewicz, B. R. Noack, F. Thiele, and G. Tadmor. 
Generalized mean-field model for flow control using continuous mode 
interpolation. In 3rd AIAA Flow Control Conference, San Francisco, 
Ca, USA, 5-8 June 2006, 2006. Invited AIAA-Paper 2006-3488. 

[40] B. R. Noack, M. Morzynski, and G. (eds.) Tadmor. Reduced- Order 
Modelling for Flow Control. Number 528 in CISM Courses and Lectures. 
Springer- Verlag, Berlin, 2011. 

[41] D. Rempfer and F. H. Fasel. Dynamics of three-dimensional coherent 
structures in a flat-plate boundary-layer. J. Fluid Mech., 275:257-283, 
1994. 

[42] J. Gerhard, M. Pastoor, R. King, B. R. Noack, A. Dillmann, M. 
Morzynski, and G. Tadmor. Model-based control of vortex shedding 
using low- dimensional Galerkin models. In 33rd AIAA Fluids Confer- 
ence and Exhibit, Orlando, Florida, USA, June 23-26, 2003, 2003. Paper 
2003-4262. 

[43] M. Bergmann, L. Cordier, and J. -P. Brancher. Optimal rotary control of 
the cylinder wake using proper orthogonal decomposition reduced order 
model. Phys. Fluids, 17:097101-1. . .21, 2005. 

[44] B. Thiria, S. Goujon-Durand, and J. E. Wesfreid. The wake of a cylinder 
performing rotary oscillations. J. Fluid Mech., 560:123-147, 2006. 

[45] M. Pastoor, L. Henning, B. R. Noack, R. King, and G. Tadmor. Feed- 
back shear layer control for bluff body drag reduction. J. Fluid Mech., 
608:161-196, 2008. 

[46] J. Weller, E. Lombardi, and A. Iollo. Robust model identification of 
actuated vortex wakes. Physica D, 238:416-427, 2009. 

[47] C. A. J. Fletcher. Computational Galerkin Methods. Springer, New 
York, 1st edition, 1984. 

[48] P. Holmes, J. L. Lumley, G. Berkooz, and C. W. Rowley. Turbulence, 
Coherent Structures, Dynamical Systems and Symmetry. Cambridge 
University Press, Cambridge, 2nd paperback edition, 2012. 



Optimal Model Identification 



39 



[49] L. Sirovich. Turbulence and the dynamics of coherent structures, Part 
I: Coherent structures. Quart. Appl. Math., XLV:561-571, 1987. 

[50] R. K. Miller and A. N. Michel. Ordinary Differential Equations. Aca- 
demic Press, 1982. 



